Stripe formation in doped Hubbard ladders 
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We investigate the formation of stripes in 7r x 6 Hubbard ladders with 4r holes doped away from 
half filling using the density-matrix renormalization group (DMRG) method. A parallelized code 
allows us to keep enough density-matrix eigenstates (up to m — 8000) and to study sufficiently large 
systems (with up to 7r — 21 rungs) to extrapolate the stripe amplitude to the limits of vanishing 
DMRG truncation error and infinitely long ladders. Our work gives strong evidence that stripes 
exist in the ground state for strong coupling (U = 12t) but that the structures found in the hole 
density at weaker coupling (U = 3t) are an artifact of the DMRG approach. 
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Two-dimensional lattice models for correlated elec- 
trons are often used to describe the properties of layered 
cuprate compounds^ Despite numerous studies there is 
an ongoing controversy about the existence of stripes in 
these systems. In a hole doped system a stripe is a do- 
main wall ordering of holes and spins. The wall is made 
of a narrow hole-rich region. The spins are antiferromag- 
netically ordered between the walls and are correlated 
with a 7T phase shift across a wall. The formation of 
stripes in the ground state has been demonstrated numer- 
ically for the t— J model on (narrow) ladders^ using the 
density-matrix renormalization group (DMRG) method. 4 
For square lattices, however, the presence of stripes re- 
mains controversial^LS because a reliable investigation 
of the ground state in the thermodynamic limit is not 
possible with the methods currently available. 

Recently, attention has turned to the two-dimensional 
Hubbard model with a local electron-electron repulsion 
U and an electron hopping term t. For U = this model 
describes a Fermi gas, which obviously has no stripes in 
the ground state. Moreover, no instability toward the 
formation of stripes has been found in the weak-coupling 
limit U -C t using renormalization group techniques £ 
In the strong-coupling limit U ^> t, however, the Hub- 
bard model can be mapped onto a t — J model with 
J = 4t 2 /U <C t, which does have stripes in the ground 
state, at least on narrow ladders with J r* 0.35tm& 
Therefore, investigating the formation of stripes in the 
Hubbard model at finite coupling U jt could significantly 
improve our understanding of these structures. More- 
over, such investigations should reveal the true capability 
of the various methods used to study stripes much better 
than calculations for the t — J model alone. 

An early DMRG investigation of 3-leg Hubbard lad- 
der^ found that stripes formed in the ground state only 
for U > 6t. In a recent DMRG calculation White and 
Scalapino^ have shown that a narrow stripe appears in 



the ground state of 6-leg Hubbard ladders (more pre- 
cisely, 7x6-site clusters) for U > 6t. For weaker couplings 
the hole and spin densities show structures which are in- 
terpreted as a broad stripe. In both works, however, no 
finite-size scaling has been performed and the amplitude 
of the hole density modulation has not been investigated 
systematically as a function of DMRG truncation errors. 
Thus it is not clear if the observed structures are really 
the signature of a striped ground state in the limit of 
infinitely long ladders; they could be finite-size effects or 
an artifact of the DMRG method. 

Here, we report on a DMRG investigation of stripes in 
6-lcg Hubbard ladders with up to 28 rungs. Keeping up 
to m = 8000 density-matrix eigenstates the amplitude of 
the hole density modulation can be extrapolated to the 
limit of vanishing DMRG truncation errors for systems 
with up to 21 rungs. This allows us to perform a reliable 
finite-size scaling analysis of the hole density modulation. 
Calculations for systems of that size are made possible 
by a parallelized shared-memory DMRG code which runs 
efficiently on current supercomputer architectures^ We 
show that the stripes found by White and Scalapinc 2 
are stable in the limit of an infinitely long ladder for 
strong coupling U = 12i. For weak coupling (U = 3t), 
however, the hole density fluctuations found in Ref.0are 
an artifact of truncation errors and boundary conditions. 

The two-dimensional Hubbard model on a R x L ladder 
is defined by the Hamilton operator 

H = —t (^x,y,cr^x,y+l,(T + ^x,y,cr^x+l,y,a + h.C.J 

x,y,cr 

+U n x , y ,\n Xj y^ , (1) 

where x = 1, . . . , R is the rung index and y = 1, . . . , L is 
the leg index, c\ a and c x ,y,a are creation and annihi- 
lation operators for an electron with spin a =|, \ at site 
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(x,y), and h x . 



st 



is the corresponding den- 



sity operator. Here we exclusively consider the Hubbard 
model on 6-leg ladders (L = 6) with R = 7r rungs for 
r = 1, . . . ,4. Cylindrical boundary conditions were used 
(closed in the rung [y] direction and open in the leg [x] 
direction), because they are the most favorable ones for 
DMRG simulations. Moreover, open boundaries break 
the translational invariance of the system, allowing spin 
and charge structures to appear as local density varia- 
tions in a finite ladder. If periodic boundary conditions 
were used, one would have to analyze correlation func- 
tions to detect stripes in finite ladders. Since we are in- 
terested in the ground state of the hole-doped regime, we 
consider a system with N = Ar holes doped in the half- 
filled band, corresponding to RL — N = 38r electrons. 
The average hole density is n = N/RL = 4/42 ks 0.095 
for all cases, as in Ref. 0- 

We employ a recently developed parallelized DMRG 
codeii to determine the ground state properties of this 
Hubbard model. Our DMRG program implements the 
standard finite-system algorithm— for two-dimensional 
lattices. Parallelization of the most time-consuming tasks 
allows us to carry out calculations of unprecedented mag- 
nitude. For the calculations presented here we have kept 
up to to = 8000 density-matrix eigenstates per block for 
systems with up to R x L = 168 sites. This requires up 
to four weeks walltime and 100 GBytes of memory per 
run on eight processors of an IBM p690 node. For com- 
parison, it takes about 6 hours to reproduce the results 
of Ref. |a for 7 x 6 clusters with m = 3600. DMRG cal- 
culations have already been performed for the Hubbard 
model on larger systems (square lattices or ladders) than 
our 28 x 6-site clusters but for a significantly smaller num- 
ber of density matrix eigenstates (to < 2000) The 
computational cost of these simulations was at least an 
order of magnitude lower than in the present work. 

For our 6-leg Hubbard ladders, the standard DMRG 
method yields the ground state energies and various ex- 
pectation values for the ground state of the system in- 
vestigated. Here, we focus on the hole density 



h{x, y) = l- (h x>y j + n x . yd ) 
and the staggered spin density 

s(x,y) = (-l) x+y (n x , y ,i - n x>yi i) 



(2) 



(3) 



where (. . .) represents the (DMRG) ground state expec- 
tation value. In the first few lattice sweeps of our DMRG 
calculations or for a small number to of density-matrix 
eigenstates per block (m < 1000), the DMRG wavefunc- 
tion reaches a 'metastable' state^^S which depends es- 
sentially on the initial conditions, i.e., on the detail of the 
method used to construct the lattice in the first sweep 
(for more detail about the DMRG method, see Ref. 3. 
The hole and staggered spin densities show irregular fluc- 
tuations in both the rung and the leg directions at that 
point of the DMRG calculation. 

For all system sizes and coupling strengths investi- 
gated, the DMRG wavefunction "tunnels" to a stable 
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FIG. 1: (Color online) DMRG ground state energy per site 
versus m for the 21 x 6 ladder at U = 12t, using (circles) 
and not using (squares) block reflection in the middle of the 
lattice. Every data point corresponds to one sweep of the 
DMRG procedure. 




FIG. 2: (Color online) Hole (circle) and staggered spin 
(square) densities in the leg [x] direction calculated on a 21 x 6 
Hubbard ladder with 12 holes for U = 3t (open symbols) and 
U = 12t (solid symbols). 



state after several sweeps and for sufficiently large m, 
as reported in Ref. for 7x6 clusters. We have observed 
that the larger the system length R, the larger must m 
be to reach the stable state, ranging e.g. from m « 600 at 
R = 7 to to w 2200 at R = 21 in the U = 12i case. This 
state is then essentially independent of the initial condi- 
tions, but it is nevertheless essential to make m as large 
as feasible in order to get sufficient data for a reliable 
extrapolation of observables (see below). The tunneling 
occurs for smaller m when it is possible to utilize the 
block reflection technique (see Fig. QJ. Note, however, 
that the combination of spin and charge density fluctu- 
ations can easily break the symmetry between left and 
right DMRG blocks and that the block reflection tech- 
nique should not be used in that case as it can lead to 
incorrect results. For the 7x6 case, it is interesting to 
note that the transition occurs at a significantly smaller 
to with our algorithm than in Ref. [2J. This may also 
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be attributed to a more favorable choice of initial condi- 
tions. Our results for large m are in full agreement with 
those presented in Ref. 0- The "tunneling" is marked by 
a sharp drop in energy, while the spin and hole densi- 
ties become more regular. In particular, the hole density 
and the staggered spin density are almost constant in 
the rung direction. The stability of this DMRG 'ground 
state' is demonstrated by the systematic behavior of the 
energy (for all system sizes) and expectation values (for 
systems with up to 21 rungs) as a function of the dis- 
carded weight (see the discussion below). Some of our 
results for 28 x 6 ladders are inconclusive because an in- 
sufficient number of sweeps has been performed for the 
charge and spin density profiles to reach convergence. 

On the 7r x 6 ladders with 4r holes investigated here, 
r stripes with 4 holes each appear in the DMRG ground 
state. These stripes are clearly seen in the hole density 
modulation in the leg direction 



(4) 



which is shown in Fig. [21 for a 21 x 6 ladder with U = 
12t. In the same figure, one sees that the staggered spin 
density in the leg direction 



(5) 



is finite and changes sign exactly where the hole den- 
sity h(x) is maximal. Therefore, the specific features of 
stripes are clearly observed in the DMRG ground state 
densities. Note, however, that the finite staggered spin 
density is an artifact of our DMRG method^, which does 
not use the full spin symmetry. In the true ground state 
of a finite ladder one expects s(x,y) = 0. In Fig. [5] the 
results for U = M appear qualitatively similar to those 
obtained for U — 12t although the amplitudes of the 
density fluctuations are smaller for the weaker coupling. 
Nevertheless, one notices that the hole and spin density 
profiles for U — 3t are less regular than for U = 12t. 

To make a quantitative analysis of these structures we 
have carried out a systematic spectral analysis of the hole 
and staggered spin densities. The spectral transforms are 
defined as 



F{k x ,k y ) = J L{R+1) E Mk x xy k yVf(x, y) (6) 



with k x = z x n/(R + 1) for integers z x = 1, . . . , R and 
k y = 2-KZyjL for integers —L/2 < z y < L/2. Here f(x, y) 
and F{k x ,k y ) represent either the hole density h{x,y) 
and its transform H(k x ,k y ) or the staggered spin den- 
sity s(x, y) and its transform S(k x , k y ). The transforma- 
tion in the rung direction (with periodic boundary condi- 
tions) is the usual Fourier transform. In the leg direction 
(with open boundary conditions) we use an expansion in 



particlc-in-thc-box eigenstates because this is a natural 
basis for a finite open system. In the infinite-ladder limit 
R — > oo this transformation becomes equivalent to the 
standard Fourier transformation. As the systems consid- 
ered have a reflection symmetry (x — > R+l — x), the hole 
spectral transform H{k x , k y ) always vanishes for even in- 
tegers z x while the spin spectral transform S(k x , k y ) van- 
ishes for odd z x if R is odd and for even z x if R is even. 
Moreover, in the converged DMRG ground state we ob- 
served uniform behavior of h(x,y) and s(x,y) along the 
rungs. This implies that the spectral weight is concen- 
trated at k y = for both densities. 

In Fig. [3] we show the power spectrum (squared norm 
of the spectral transforms) of the hole and staggered spin 
densities presented in Fig. [21 In both cases, the power 
spectrum has been normalized by its total weight 



\F(k x , ky 



(7) 



which we denote H 2 and S 2 for the hole and spin power 
spectrum, respectively. For U ~ \2t one sees that both 
hole and spin power spectra have a single strong peak 
containing most of the spectral weight (92 % and 84%, 
respectively). For U = 3t, however, we observe a broad 
distribution without a clearly dominant mode k x in the 
hole power spectrum. We find similar results for other 
ladder lengths R < 21. For R = 28 the power spectra 
are mostly inconclusive because of the non-convergence 
of the hole and spin densities mentioned previously. 

For U = 12t the observed positions of the dominant 
peaks in the hole and spin spectral transforms for r <E 
{1, 2, 3} can be extrapolated to the R = oo limit, yielding 
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respectively, which agrees perfectly with the expected 
values corresponding to one stripe every seven rungs in 
an infinitely long ladder. For U = 3i, the position k x 
of the maximum in the spectral transforms is not always 
well defined (for instance, it changes with the number 
m of density matrix eigenstates kept even for large m) 
and thus a quantitative analysis of k x for R — * oo is not 
possible. 

All DMRG calculations suffer from truncation errors 
which are reduced by increasing the number m of retained 
density matrix eigenstates (for more details, see Refs.0). 
The error in the ground state energy is proportional to 
the discarded weight W m which is defined as the total 
weight of the discarded density-matrix eigenstates: 
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2=771+1 



(10) 
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FIG. 3: (Color online) Normalized power spectrum of the hole 
(solid bar) and staggered spin (open bar) densities presented 
in Fig. for a 21 x 6 Hubbard ladder with U = 3t (upper 
panel) and U = 12t (lower panel). 



FIG. 4: (Color online) Maximum if max of \H(k x , 0)| in a 21 x 6 
system with U = 12t (solid circle) and U = 3t (open circle) as 
a function of the square root of the discarded weight W m for 
various numbers of density-matrix eigenstates 1800 < m < 
8000. Dashed lines are linear fits. 



Here, d is the dimension of the density matrix and Wi 
is its ith eigenvalue. Thus one can achieve a greater ac- 
curacy and obtain an estimate of the error with a linear 
extrapolation of the DMRG energy to the limit of van- 
ishing truncation errors W m — > (for an example of this 
extrapolation, see Ref . ITo|) . Using this technique we have 
found that the error in the ground state energy per site 
is typically about 4 x 10~ 3 £ (U = 12t) and 7 x lCT 3 i 
(U = 3t) for the largest number of density-matrix eigen- 
states kept (to = 8000 and to = 6000, respectively) and 
R < 28. Consequently, the error in the total energy is 
of the order of t for the largest system (168 sites). The 
energy separation between ground state and the lowest 
excited states, which is of the order of a small fraction 
of t, is thus significantly smaller than the error in the 
total energy. Therefore, the DMRG wavefunctions ob- 
tained in our calculations are not accurate descriptions of 
the true ground states. Although the DMRG wavefunc- 
tions converge systematically to the true ground states 
(as shown by the linear scaling of the energy with W m ), 
for to < 8000 they still have significant overlaps with 
other eigenstates. Expectation values calculated with 
such a DMRG wavefunction (i.e., for a given to) could 
thus be quite inaccurate. In order to get reliable results 
we will in the following carefully analyze the scaling of 
observables with increasing to. 

If a variational ground-state wavefunction, as used in 
the DMRG algorithm, is known up to an error of s, the 
corresponding error in the energy is of the order of e 2 . 
Other observables, whose operators are nondiagonal in 
the energy basis, converge only with e. For DMRG we 
know that e 2 cx W m (cf. Ref.llOr). thus expectation values 
of operators are polynomials of \JW m for small discarded 
weights W m . For the maximum -ff max of the absolute 
hole spectral transform \H(k x ,0)\ we find a linear scal- 
ing with \fW m (see Fig. . Such a scaling has already 
been found for other density modulation amplitudes. 13 



This allows us to extrapolate -ff max to the limit of van- 
ishing truncation errors. We note that ff max decreases 
with decreasing W m . This corresponds to a diminution 
of the stripe amplitude with increasing to (i.e., decreas- 
ing W m ) for sufficiently large to > 2000. This diminution 
can be seen directly in the hole density profile h(x) (for 
instance, in Fig. 3b of Ref. 0). For U = 12t the ex- 
trapolated values of -ff max are clearly finite as shown in 
Fig. 0] for a 21 x 6 ladder. Thus we conclude that the 
hole density fluctuations found on finite ladders are not 
an artifact of DMRG truncation errors but a feature of 
the true ground state for U = 12t. For U = 3t, ff max 
extrapolates to very small values as W m — > 0. Typically, 
the range of yJW m over which we observe a linear behav- 
ior in y/Wm is smaller than the smallest value of \JW m 
used in the extrapolation. This can be seen for the exam- 
ple shown in Fig. 01 The uncertainty in the extrapolated 
H ma x is thus larger than its value for U = it. There- 
fore, the hole density fluctuations could be the result of 
DMRG truncation errors and the true ground state could 
be uniform in that case, i.e., i? max = for W m — > 0. 

Extrapolating the maximum 5 max of the absolute spin 
spectral transform 15(^,0)1 to the limit W m — > turns 
out to be more difficult than extrapolating H max . Con- 
trary to -ff max , S max has not reached an asymptotic 
regime (as a function of W m ) for the largest number to of 
density matrix eigenstates we have used. This difference 
is probably due to the smaller energy scale of spin excita- 
tions compared to that of charge excitations, resulting in 
a DMRG ground state which describes the charge prop- 
erties of the true ground state correctly but not its spin 
properties. Nevertheless, for the smallest (7 x 6) ladder 
an extrapolation of S max to W m — ► using linear and 
quadratic fits suggests that S max vanishes for W m — ► 
and thus that the true ground state has no spin density 
fluctuations, s(x,y) = 0, as expected (sec Fig. [5}. The 
artificial breaking of the spin symmetry is similar for all 
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FIG. 5: (Color online) Maximum Smax of |5(fe,0)| for the 
7x6 system. Dashed line: linear fit, dotted line: quadratic 
fit, solid line: quadratic fit with constraint S max (0) = 0. 
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FIG. 6: (Color online) Amplitude H max /*/R~ oc ho of the hole 
density modulation for a fixed number (6000 < m < 8000) of 
density-matrix eigenstates (square) and extrapolated to the 
limit Wm — » (circle) as a function of the inverse ladder 
length 1/R for U = 12t (solid symbols) and U — 3t (open 
symbols). Dashed lines are linear fits. 



couplings U, indicating that it does not affect the forma- 
tion of stripes (i.e., the existence of a hole density modu- 
lation) , which is a strongly U -dependent phenomenon as 
shown here. 

Typically, the discarded weight W m is about 10~ 5 or 
smaller for m = 8000. Although this appears to be a 
small value, the above discussion shows that errors in 
the ground state energy and H max are still quite large for 
that number m. This confirms that the absolute value of 
the discarded weight W m alone does not give a reliable 
estimate for errors on physical quantities. 



A ladder with a periodic array of stripes has a modu- 
lation of the hole density (charge density wave) 



h(x) 

which corresponds to 



ho sm(k^x) 



\H(k?,0)\ 



(R + l)L 
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If the amplitude ho of the hole density modulation is fi- 
nite in the limit of an infinitely long ladder (R — > oo), the 
maximum H max of the spectral transform must diverge 
as y/M for R — * oo. In Fig.Elwe show the finite-size scal- 
ing of -ffmax/ VR ~ ho as a function of the inverse ladder 
length for U = 3t and U = 12t. The values of H max /\^R 
obtained for a fixed number m of density matrix eigen- 
states converge to finite values for R — > oo, suggesting 
the existence of stripes in this limit for both couplings 
U. After extrapolation to the limit of vanishing DMRG 
truncation errors, however, H max /\/R seems to vanish 
for large R in the case U = 3i while it still converges to a 
finite value for U = 12t (see Fig. [SJ|. This indicates that 
stripes exist in the ground state of infinitely long lad- 
ders for sufficiently strong coupling such as U — I2t but 
that the hole and spin structures found in finite ladders 
for weak couplings such as U = 3t are artifacts of open 
boundaries and DMRG truncation errors. It has already 
been observed in other problcmsi^ that DMRG trunca- 
tion errors can result in an artificial broken symmetry 
ground state after extrapolation to infinite system sizes, 
while extrapolating first to the limit of vanishing trun- 
cation errors restores the true ground state symmetry in 
the thermodynamic limit. 

In conclusion, we have investigated the formation of 
stripes in 6-leg Hubbard ladders at a hole doping of 9.5%. 
Using a parallelized DMRG code we have been able to 
determine the amplitude of the hole density modulation 
in the limits of vanishing DMRG truncation errors and 
infinitely long ladders. Our results show that stripes exist 
in the ground state of these systems for strong but not 
for weak couplings. 
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